install.packages('https://cran.r-project.org/src/contrib/Archive/interflex/interflex_1.1.3.tar.gz', 
                 repos=NULL, 
                 type='source')


library(magrittr)
library(stringr)
library(plyr)
library(tidyverse)
library(interflex)
library(gridExtra)
library(readxl)


rm(list=ls())
home = 'C:/Users/jdt34/Dropbox/VNA_Responsiveness/Short Article/JOP-dataverse/'


permutation.results = paste0(home, 'Data/RI-analyses.Rds') %>%
  readRDS
experimental.results = paste0(home, 'Data/experimental-analyses.Rds') %>%
  readRDS


dv_survey = paste0(home, 'Data/survey-outcomes.xlsx') %>%
  read_xlsx %>%
  mutate(Treatment=factor(x=Treatment, 
                          levels=c('Control',
                                   'Citizen',
                                   'Firm')),
         Missing=as.integer(is.na(Q1)),
         CitizenXProp.Citizen=Citizen*Prop.Citizen, 
         FirmXProp.Firm=Firm*Prop.Firm) %>%
  subset(!is.na(Treatment))
class(dv_survey) = 'data.frame'
colnames(dv_survey)[colnames(dv_survey)=='Q1'] = 'Outcome'


dv_pooled = paste0(home, 'Data/pooled-outcomes.xlsx') %>%
  read_xlsx %>%
  mutate(Treatment=factor(x=Treatment, 
                          levels=c('Control',
                                   'Citizen',
                                   'Firm')),
         CitizenXProp.Citizen=Citizen*Prop.Citizen, 
         FirmXProp.Firm=Firm*Prop.Firm) %>%
  subset(!is.na(Treatment))
class(dv_pooled) = 'data.frame'
colnames(dv_pooled)[colnames(dv_pooled)=='Spoke'] = 'Outcome'


tmp = subset(dv_survey, !is.na(Outcome)) %>%
  ggplot(aes(x=Prop.Citizen, fill=factor(Citizen))) +
  geom_histogram(color='black', bins=50) +
  scale_fill_manual(values=c('white','black'))
citizen.survey.histogram = layer_data(tmp)[,c('fill','xmin','xmax','ymin','ymax')] %>% 
  mutate(xmax=0.974*(xmax+abs(min(xmin))), xmin=0.974*(xmin+abs(min(xmin))), 
         ymin=ymin/160-1, ymax=ymax/160-1)


tmp = ggplot(data=dv_pooled, 
             aes(x=Prop.Citizen, fill=factor(Citizen))) +
  geom_histogram(color='black', bins=50) +
  scale_fill_manual(values=c('white','black'))
citizen.pooled.histogram = layer_data(tmp)[,c('fill','xmin','xmax','ymin','ymax')] %>% 
  mutate(xmin=0.974*(xmin+abs(min(xmin))), 
         xmax=0.974*(xmax+abs(min(xmin))), 
         ymin=ymin/160-1, 
         ymax=ymax/160-1)
rm(tmp)


rimer = function(result, X, Z, dv, Label) {
  me.sigs = NULL
  Zrange = seq(from=0, 
               to=0.71, 
               by=0.01)
  if(result=='tab2col3') {
    observed_data = dv_survey
    TreatedHistogram = subset(citizen.survey.histogram, 
                              fill=='black')
    ControlHistogram = subset(citizen.survey.histogram, 
                              fill=='white')
  } else if(result=='tab3col3') {
    observed_data = dv_pooled
    TreatedHistogram = subset(citizen.pooled.histogram, 
                              fill=='black')
    ControlHistogram = subset(citizen.pooled.histogram, 
                              fill=='white')
  }
  ri.ests = (subset(permutation.results, 
                    Result==result & term==X)$estimate + 
               (subset(permutation.results, 
                       Result==result & term==paste(X, Z, sep=':'))$estimate %*% t(Zrange)))
  our.est = (subset(experimental.results, 
                    Result==result & term==X)$estimate + 
               subset(experimental.results, 
                      Result==result & term==paste(X, Z, sep=':'))$estimate*Zrange)
  for(z in 1:length(Zrange)) {
    qest = ecdf(ri.ests[,z])
    me.sigs[z] = qest(our.est[z]) %>%
      multiply_by(1e3) %>%
      round %>%
      as.integer %>%
      mapvalues(from=seq(0L, 1000L, 1L),
                to=c(rep('p \u2264 0.01', 6), 
                     rep('p \u2264 0.05', 20), 
                     rep('p > 0.05', 950), 
                     rep('p \u2264 0.05', 20), 
                     rep('p \u2264 0.01', 5)),
                warn_missing=F)
  }
  me.sigs = factor(x=me.sigs, 
                   levels=c('p > 0.05',
                            'p \u2264 0.05',
                            'p \u2264 0.01'))
  Labels = list(expression(p > 0.05), 
                expression(p <= 0.05), 
                expression(p <= 0.01))
  if.D = X
  if.X = Z
  if.Z = c('FirmXProp.Firm', 
           setdiff(subset(experimental.results, 
                          Result==result)$term, 
                   c(if.D, if.X, '(Intercept)','Citizen:Prop.Citizen','Firm:Prop.Firm')))
  kern.est = inter.kernel(Y='Outcome', 
                          D=if.D, 
                          X=if.X, 
                          Z=if.Z, 
                          data=observed_data, 
                          na.rm=T, 
                          CI=F, 
                          cores=6)$est[[1]][,c('X','ME')] %>%
    set_colnames(c('Z','dYdX'))
  ri.ests = as.data.frame(ri.ests) %>%
    mutate(ID=1:1e4) %>%
    gather(key='Z', value='dYdX', -ID) %>% 
    mutate(Z=mapvalues(x=Z, 
                       from=paste0('V', 1:length(Zrange)), 
                       to=Zrange),
           Z=as.numeric(Z))
  data.frame(Significance=me.sigs, 
             Z=Zrange, 
             dYdX=our.est) %>%
    ggplot(aes(x=Z, y=dYdX, fill=Significance)) +
    geom_line(data=ri.ests, 
              aes(x=Z, y=dYdX, group=ID), 
              color='black', alpha=0.01, inherit.aes=F) +
    geom_line(data=kern.est, 
              aes(x=Z, y=dYdX), 
              color='#e63e40', linetype=1, alpha=1, size=0.5, inherit.aes=F) + 
    geom_hline(yintercept=0, linetype=2, color='white') +
    geom_point(size=2, shape=21) +
    geom_rect(data=TreatedHistogram, 
              aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
              fill='black', color='black', size=0.25, inherit.aes=F) +
    geom_rect(data=ControlHistogram, 
              aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax), 
              fill='white', color='black', size=0.25, inherit.aes=F) +
    geom_text(data=data.frame(x=0.05, y=0.875, label=Label, stringsAsFactors=F),
              aes(x=x, y=y, label=label),
              inherit.aes=F, hjust=0, vjust=0.5, size=5) +
    scale_x_continuous(breaks=seq(from=0, 
                                  to=1, 
                                  by=0.2), 
                       labels=paste0(str_pad(string=seq(from=0, 
                                                        to=100, 
                                                        by=20), 
                                             width=3, 
                                             side='left', 
                                             pad=' '), 
                                     '%')) +
    scale_fill_brewer(labels=Labels, type='seq', palette='YlGnBu', direction=1, drop=F) +
    guides(fill=guide_legend(override.aes=list(size=2))) +
    labs(x='Moderator: % citizen', 
         y='Marginal effect of citizen treatment') +
    coord_fixed(ratio=0.25, xlim=c(0,1), ylim=c(-1,1), expand=F) +
    theme_minimal() +
    theme(axis.text=element_text(size=14, color='black'),
          axis.title.x=element_text(size=14, color='black'),
          axis.title.y=element_text(size=14, color='black', hjust=1),
          strip.text=element_text(size=9),
          legend.title=element_text(size=14),
          legend.text=element_text(size=14),
          legend.position='bottom',
          plot.margin=unit(c(0.1, 0.3, 0, 0.1), 'in'),
          panel.border=element_rect(color='black', fill=NA))   
}


fig_1_left = rimer(result='tab2col3', 
                   X='Citizen', 
                   Z='Prop.Citizen', 
                   Label='A) Prepared for debate')
fig_1_right = rimer(result='tab3col3', 
                    X='Citizen', 
                    Z='Prop.Citizen',
                    Label='B) Spoke in any forum')
fig_1 = grid.arrange(fig_1_left, fig_1_right, nrow=1)


ggsave(filename='figure-1.png', 
       plot=fig_1, 
       path=paste0(home, 'Figures/'), 
       width=12, 
       height=3.75, 
       units='in')
ggsave(filename='figure-1.eps', 
       plot=fig_1, 
       path=paste0(home, 'Figures/'), 
       width=12, 
       height=3.75, 
       units='in', 
       device=cairo_ps)
